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ABSTRACT 

This review covers the theory and application of spectral collocation 
methods. Section 1 describes the fundamentals, and summarizes results per- 
taining to spectral approximations of functions. Some stability and con- 
vergence results are presented for simple elliptic, parabolic and hyperbolic 
equations. Applications of these methods to fluid dynamics problems are 
discussed in Section 2. 
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INTRODUCTION 


Spectral collocation methods form an efficient and highly accurate class 
of techniques for the solution of nonlinear partial differential equations* 
They are characterized by the expansion of the solution in terms of global 
basis functions, where the expansion coefficients are computed so that the 
differential equation is satisfied exactly at a set of so-called collocation 
points* The fundamental unknowns are the solution values at these points* 
The expansion is used only for the purpose of approximating derivatives* 

The popularity of these methods arises from several advantages which they 
have over common finite difference methods* First, they have the potential 
for rapidly convergent approximations. The dual representation in physical 
and transform space allows for automatic monitoring of the spectrum of a solu- 
tion, providing thus a check on resolution. If certain symmetries exist in 
the solution, spectral methods allow the exploitation of these symmetries* 
Finally, the methods have low or non-existent phase and dissipation errors. 
For these reasons, spectral methods have become the primary solution technique 
in such areas of computational fluid dynamics as the simulation of turbulent 
flows and the computation of transition to turbulence. They are becoming 
increasingly viable in certain areas of computational aerodynamics such as 
compressible flows, boundary layers and heat transfer* 


1.1 BASICS 

The spectral representation of a function or solution to a differential 
equation is an expansion in a series of global orthogonal polynomials: 
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u(x) = T a U (x). (1) 

n n 

n=-oo 

If u is periodic, the usual basis (or expansion) functions, U^, are com- 
plex exponentials, e inx . If u is not periodic, Chebyshev or Legendre poly- 
nomials are often used. In general, expansions in terras of polynomials which 
are solutions of singular Sturm-Liouville problems are appropriate. This is 
discussed with regard to spectral solutions of a variety of boundary value 
problems in Gottlieb and Orszag [1] and in a more formal setting by Quarteroni 
[2]. Here we summarize the essential results for the one-dimensional case; 
the basis functions for multi-dimensions are defined by tensor products of 
one-dimensional basis functions. 

The formally self-adjoint Sturm-Liouville eigenvalue equation is written 
as 

LU = - {(plT)' + qU } = X U (2) 

n w 1 *n M n J nn 

x€[a,b] 

where we assume that p(x) >0 is differentiable and q(x) ^ 0 is contin- 
uous. The function w(x) >0 is the weight function. It is well known 
(e.g., Courant and Hilbert [3]) that with the proper boundary conditions the 
eigenvalues of this equation form a countably infinite sequence. Also, the 

set of eigenfunctions {U n } is orthogonal and forms a basis for the Hilbert 
2 

space L w with the weighted inner product 

b 

(u,v) = / uvwdx (3) 

w J 
a 

1/2 

and norm (ul * (u,u) . Since the U_, are orthogonal, the coeffi- 

W W 11 
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cients a n of equation (1) are defined by the relation a n = ( u , U n ) w / H H w . 

By choosing v = LU n in equation (3), u = U n and integrating by parts, 

a = -J— (Lu,U ) + -J- [P(lTu - U u')]J* . (4) 

n A n w X n n a 

n n 


If the boundary conditions are periodic or if u and the U n exactly 
satisfy homogeneous boundary conditions, then the second term is zero and 

I 

l 


n 


n 


(5) 


From this point, Quarteroni [2] iterates to derive 


a = — (u, x,U ) 
n m (m) n w 
A 

n 


( 6 ) 


where U( m ) = Lu^.^. From the properties of the X n , as n+». 


a < — t — I u, vl . 
1 n 1 — 2m (m) w 
n 


(7) 


Thus, the magnitude of the coefficients of the expansion (1) depends only on 
the degree of smoothness of the solution. In particular, if u and the co- 
efficients p and q are infinitely differentiable, then the a n decay 
faster than any polynomial of 1/n. 

To ensure that inhomogeneous boundary conditions do not affect the rate 
of convergence for non-periodic problems, the second term of equation (4) must 
always be zero. This occurs if p(x) is identically zero at the bounda- 
ries. The equation (2) is called a singular Sturm-Liouville equation if p(x) 
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= 0 at least at one of the boundaries. In general, one needs p(x) =0 at 
both boundaries. The polynomials for which this is true are known as Jacobi 
polynomials. These polynomials are orthogonal with respect to the weight 
w( x ) = a(l - x) a (l + x) V and p(x) = 6(1 - x) 0+1 (l + x) v+1 where 
a, 3, a, v are real constants. The most commonly used polynomials for 
spectral computations are the Chebyshev polynomials where a = v = -1/2 
and a = 3 = 1 and the Legendre polynomials where o = v = 0 and a = 3 = 1. 

In practice, the choice of basis depends on the type of boundary condi- 
tions and the efficiency with which the expansions can be computed. Fourier 
methods are certainly appropriate for problems with periodic boundary condi- 
tions. They can also be computed efficiently by using fast Fourier transform 
techniques (see Gottlieb and Orszag [1]). For non-periodic cases, Chebyshev 
spectral methods are popular because they can use a fast cosine transform. 
Legendre methods are not as commonly used because they lack a fast transform 
method. 


1.2 COLLOCATION METHODS 

The application of spectral methods consists of approximating the infi- 
nite sum in equation (1) with a finite sum by projecting the function onto a 

p 

finite subspace of L^. In particular, the collocation projection consists 
of finding the polynomial which takes on the exact value of the original func- 
tion at a finite number of grid (or collocation) points. In other words, the 
collocation projection is interpolation. This method, often termed pseudo- 
spectral, is particularly useful for nonlinear problems. 
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| 

1.2.1 FOURIER COLLOCATION 

For problems with periodic boundary conditions, the Fourier expansion of 
j a function u(x) is given by the infinite series 

| 00 

u(x) = l a e ikx . (8) 

k=-» K 


The collocation projection is defined by the discrete Fourier transform pair 

I 

I 


V 


N/2-1 

I 

k=-N/2 


V 


ikx 


(9) 


where the coefficients are defined by 



N-l -ikx 

I u ( x . )e 
j=0 3 



2 


+ 1 


> 


N 

2 


1 . 


The collocation points, x j , are uniform on the interval [0,2ir] 


( 10 ) 


Xj = 2irj/N j = 0, 1 ,2, . . .N-l . 


(ID 


The transforms (9) and (10) are almost always computed by the use of a fast 
Fourier transform if N is a highly composite integer such as N = 21*39. 

Derivatives of the function u at the collocation points are approxi- 
mated by the derivatives of the interpolating polynomial. Thus, the £,th 
derivative of u is approximated 


d *V 


N/ 2-‘ ikx 

l (ik) u, e 
k=-N/2 


( 12 ) 
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From the form of equation (12), it is clear that the evaluation of the deriva- 
tive at the collocation points can also be computed efficiently with a fast 
Fourier transform. See Hussaini, Streett, and Zang [4] for more information 
on the implementation of the Fourier collocation method. 

1.2.2 CHEBYSHEV COLLOCATION 

The collocation points for using a Chebyshev method to approximate a non- 
periodic function are usually defined by 

Xj = -cos(irj/N) j = 0,1, ...,N. (13) 

These points are the extrema of the N 1 -* 1 order Chebyshev polynomial, T N (x), and 
are obtained from the Gauss-Lobatto integration formula (see Davis and 
Rabbinowitz [5]). 

The collocation projection operator is defined as the interpolation 

N 

P XT u = J a T (x) (14) 

N **_ n n 
n=0 

where the coefficients are defined by 

(2 j = 1 ,N 

where C . = \ • (15) 

( 1 otherwise 

Again, derivatives of u at the collocation points are approximated by 
the derivative of the interpolating polynomial evaluated at the collocation 
points. The first derivative, for example, is defined by 


2 1 ? u ( x )T (x ) 

n J 3 
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~ = y T (x) 

dx 4 a n V x; 

n=0 


■Si! - <* 


(16) 


a (1) = o 

a N U 


C a^ = a^ + 2(n + l)a . 
n n n+2 n+1 


a N 1)1 |**«jO • 


The transform pair given by equations (14) or (16) and (15) can be 
efficiently computed with a fast cosine transform. Equivalently, the inter- 
polating polynomial and its derivatives can be computed using matrix multipli- 
cation. The matrices for the Chebyshev collocation method are conveniently 
collected in the review by Gottlieb, Hussaini, and Orszag [6]. For N < 32, 
this approach is competitive with using a fast cosine transform, at least on 
serial computers. 


1.3 APPROXIMATION THEORY (COLLOCATION) 

1.3.1 FOURIER COLLOCATION 

The problem of how well P N u approximates u for Fourier approximations 
has been discussed by Kreiss and Oliger [7], Pasciak [8], and by Canuto and 
Quarteroni [9]. See also Mercier [10]. It is most convenient to express the 
interpolation results in terms of a Sobolev space, H m (0,2ir). This is a 
Hilbert space with the norm 



nun 


(17) 


q 


j=o 


2 

i 


defined in terms of the seminorms 


2 

P 


l 


k=-°° 



(18) 


The use of the discrete Fourier transform pair (9), (10) represents the 
projection of the Sobolev space onto the space S n (0,2tt), the space of 
Fourier polynomials of degree N. 

The primary interpolation result is given by Theorem 1: 


Theorem 1: For any 0 £ p _< q with q > 1/2 there exists a 

constant C independent of u and N such that 


Hu - P„ull < CN p_q |u| . 
N p — 1 1 q 


(19) 


Proof : See Pasciak [8]. 

Equation (19) states that the rate of convergence depends (through the 
order of Sobolev norms) only on the smoothness of the function being approxi- 
mated. This type of error decay is known as spectral accuracy. In practice, 
one sees errors which decay exponentially and hence spectral accuracy is often 
called exponential accuracy. Several applications described in Section 2 
exhibit exponential accuracy. The term infinite order accuracy is also used 
often to refer to the case as q ». 
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Exponential accuracy has been shown explicitly by Tadraor [11] for func- 
tions u which are also analytic in the complex plane. 

t 

Theorem 2: Let u(x) be_ 2tt - periodic and analytic in a strip of 

i 

width 2sq. Then for any 0 < s < Sq 

II u - P N ull p _< CM(s)/sinh(s)N p e“ Ns (20) 

! 

where C depends on p and 

M(s) = Max |u(z) | . 

| Iraz | <_ s 

Proof ; See Tadmor [11]. 

If the solution is not very smooth, then the approximation may not be 
very good. In fact, if the function is discontinuous, the interpolant shows 
global oscillations (Gibbs phenomenon) and the approximation error decay is 
globally only first order. Smoothness is not usually a problem with the solu- 
tions of many elliptic or parabolic equations, but discontinuities are 
characteristic of the solutions of hyperbolic equations. 

It is still possible to obtain spectrally accurate approximations to non- 
smooth functions, at least away from any discontinuities, but some type of 
filtering is required. Two papers which address this issue are Majda, 
MacDonough, and Osher [12] and Gottlieb and Tadmor [13]. The first approach 
used to smooth discontinuous solutions was that of Majda, MacDonough and Osher 
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[12] whose results show that spectral accuracy can be retained if Fourier 
space filtering is applied. Since the main results refer directly to the 
solutions of hyperbolic partial differential equations, they will be discussed 
in the next subsection. 

Gottlieb and Tadmor [13] have taken the approach of smoothing in real 
space to allow the accuracy to depend on the local smoothness of the func- 
tion. The smoothing procedure consists of convoluting the collocation approx- 
imation with a regularization kernel which is localized in space. If we 
call P u the smoothed approximation to the originally oscillatory inter- 
polant P N u, the convolution takes on the form 


Pu(x) 



N-l ft _ 

l P N u(y i )’T ’ P U - y.) 
j-0 w 3 J 


( 21 ) 


where 


-•Ke- (, a < eT ) 


sin((p + •-)£) 
sin(y/0) 


( 22 ) 


is the Dirichlet kernel localized in space by the cutoff function 


e a? 2 /(5 2 -l) 

0 


Id < i 

otherwise 


(23) 


The function p ensures that the kernel does not interact with any regions 
of discontinuity. For example, for a single discontinuity at x = ir, they 
choose 0 = it|x - it [ . With this smoothing, they show that the error depends 
only on the smoothness of the cutoff function p(C): 
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Theorem: Let p (5 ) be a C^ s cutoff function satisfying 

p(0) = 1 and having support in [ — it ,tt ] . Then for any x _iii [0,2ir] 

the smoothed function Pu satisfies the estimate 

|Pu(x) - u(x) | < C II p I Max|D k u(y)|(l + 0 -2s )N 1-s s > 1. (24) 

” s |y-x| <9 tt 

0 < k £ 2s 

Proof: See Gottlieb and Tadmor [13]. 


1.3.2 CHEBYSHEV COLLOCATION 

To study the approximation properties of the Chebyshev projection (14), 

it is practical to work in a weighted Sobolev space with weight w(x) = ( 1 - 

x2)“l/2. Defining the weighted L 2 norm by 

w 

1 

II ull Z n = (u,u) = / u 2 wdx (25) 

U y W W m 


and the Sobolev norm by 


luH 


,d i u„ 2 


/ H T 1 ' n 

i-l dx 1 °’ w 


(26) 


the spectral approximation result is given by 


Theorem 4: Let q >1/2 and 0 < p £ q. 

constant C such that for all u in H^(-l,l) 

■ ■ - W 


Then there exists a 
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II u - P„ull < CN 
N p,w — 


2p-q 


lull 


q,w 


(27) 


Proof : See Canuto and Quarteroni [9]. 

So, like the Fourier approximation, the Chebyshev interpolation gives 
spectral accuracy; that is, the accuracy depends only on the smoothness of the 
function to be interpolated* Exponential convergence has also been proved by 
Tadraor [11]* This time, the function u must be analytic in an ellipse with 
foci at -1 and 1: 

Theorem 5: Assume u(x) is analytic in [-1,1] and has a regularity 

ellipse whose sum of its semi-axes equals r^ = expCrig) > 1. Then for any 
n ,0 < x] < n q we have 

1/2 

llu(x) - P ull < 8M(n )(-- y t - h(N ^ ) Ne -Nn (28) 

H* e P - 1 

where the norm is defined by 


Hull = l (l+p) 2 |u | 2 . 
H* p=0 P 


(29) 


Proof: See Tadraor [11]. 


If the function which is being approximated is discontinuous, it is still 
theoretically possible to recover a spectrally accurate solution [13] by fil- 
tering in physical space. The procedure is the same as the smoothing proce- 
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dure for the Fourier case, but the Dirichlet kernel is replaced by 


K (y) - 7 A 

p pit y 


T ;(y) 


( 30 ) 


1.4 THEORY OF SPECTRAL COLLOCATION METHODS FOR PDE'S 

Proofs of the convergence of spectral approximations to partial differen- 
tial equations are usually accomplished using energy methods which mimic 
proofs of the well-posedness of the original equations. Consequently, it is 
most convenient to discuss stability and convergence with respect to the three 
major types of partial differential equations separately. 


1.4.1 ELLIPTIC EQUATIONS 

Theoretical analysis of the convergence of Fourier collocation methods is 
simplified because of periodic boundary conditions. The elliptic problem is 
to find the function u(x) which satisfies 

Lu * f x€[0,2ir] (31) 

u(0)=u(2ir ) 


where L has the property 


(Lu,u) > a Hull 


2 

1 


a > 0. 


(32) 


The Fourier collocation approximation is obtained as described in section 
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1.2.1 and satisfies the same inequality, i.e., ¥u£Sjj 

(L u,u) > oc II u II ^ . 
c ’ — 1 

Then we have the 

Theorem 6: For t > 1 there exists a constant C such that if 

T — . O T 

f€Hp (0,2it) and u€Hp(0,2Tr) then the following estimate is optimal ; 

II u - P..UII, < C(1 + N 2 ) ( 1 _t)/ 2 ||uII . (33) 

NX— T 

Proof : See Mercier [10]. 


Chebyshev methods with both Dirichlet and Neumann boundary conditions 
have been analyzed for the elliptic differential equation of the form 


Lu = -(au ) 4* (bu) • 

xx x 


(34) 


The Chebyshev spectral collocation approximation is formally written as 


r d /n , d A dP N (bu c^ 

L c U c = " d7 (P N (a cHT U c )) + 3x + YU c 


(35) 


For Dirichlet problems, the equation is collocated at the interior points and \ 

i 

boundary conditions of the form 


u(~l) = u^ and u(+l) = u^ 


(36) 
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are specified directly at the boundary points. Stability and convergence were 
proved by Canuto and Quarteroni [14] using a variational approach. They show 


7 } bet u c be the solution to L c u c = f c where L c is defined 
by equation (35) with homogeneous boundary conditions, u r = u^ = 0 then with 
suitable conditions on a, b, a the following estimate holds 


lu - u 


i < C.N 
1 ,w — 1 


1-r, 


ull 


r,w 


+ C_N~ S || f II 


s,w 


(37) 


Proof : See Canuto and Quarteroni [14], Theorem 2.4. 

Convergence proofs for Neumann or mixed-type boundary conditions are 
available for boundary conditions applied in one of two different ways. A 
discussion of these approaches can be found in Canuto [15], [16]. The first 
approach is explicit. At interior points, the equation is collocated normally 
as in equation (35). At the boundary points, however, the collocation approx- 
imation to the derivative is written in matrix form and the boundary condi- 
tions are used to determine the value at the boundary point. Thus, the 
approximation to the boundary condition 

B,u = $u + au 
1 x 

(38) 


B u = 6u + yu 
r x 


is found by solving the system 


-16- 


N-l 

( “ + M 00 )u 0 + Sd ON U N ■ ‘ 6 I Vj 

J (39) 

+ Sd N» )U N + d N0 U 0 " Vl 

where Uj = u c (xj) and [ di j ] is the matrix for the derivative at the 

collocation points (see Gottlieb, Hussaini, and Orszag [6]). 

The convergence is very rapid for smooth solutions: 

Theorem 8: Let o > 1/2 and let u and u c be solutions to Lu = 

f and L c u c = f c where L and L c are defined as above. Then with 

explicitly applied Neumann boundary conditions the following convergence 

estimate holds 


u - u II- < CN -0 { II ull + Ilf II } 

c 2,n — 1 a+2,w a,w J 


where n * (1 ~ x^)w(x) and C is independent of N. 


Proof : See Canuto and Quarteroni [14], Theorem 3.2, 


(40) 


Canuto [16] also describes how to impose Neumann boundary conditions 
implicitly for elliptic problems. In this way, the boundary conditions are 
not exactly satisfied because what is actually solved is the modification of 
the interior approximation. For the spectral case of a pure Neumann problem, 
the first derivatives are computed normally as in equation (16). At the 
boundary points, the derivatives are replaced by the Neumann conditions. Then 
the second spectral derivatives are computed by using (16) again on the modi- 
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fied set of derivatives. This has the advantage that all of the points are 
treated the same, but the boundary conditions are not exactly satisfied. The 
boundary error does decay spectrally, however. 

theorem 9: Let u £ be the solution to L c u c = f with implicit Neumann 

boundary conditions. If u€H™(-l,l) with m > 5/2 then 

\-r-2- (±1)| < CN 4 “ m Uu8 (41) 

'9x 1 — m,w 

where C > 0 is independent of N. 

Proof: See Canuto [16]. 


The convergence in the interior is also spectral, and the estimate bounds 
both the solution and the collocation derivative. 


theorem 10: Under the assumptions of Theorem 9, 


N-l 

( 

j 


,u ‘ Vo,» + ( .i, Iu x ' <u c , x 1( *j > “j ) < 


111 ' CN 2-m iun 


m f w 


( 42 ) 


where the Wj are the Gauss-Lobatto weights at the points Xj . 


Proof: See Canuto [16] 
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1.4.2 PARABOLIC EQUATIONS 

The convergence and stability theory for linear parabolic equations, like 
the theory for elliptic equations, is fairly well developed. In particular, 
the theory has centered on studies of semi-discrete equations in which the 
spatial variation is discretized, but the time variation is left continuous. 
The emphasis, however, has been on application to boundary value problems — 
that is, on the convergence of Chebyshev collocation methods. In this section 
we survey theoretical results for initial-value problems of the form 


u„ = (Au ) + Bu + Cu + f 

t xx x 

(43) 


u(x,0) « u Q (x) 


where A, B, and C are n x n matries. The general collocation approxima- 
tion to the first, third, and fourth terms of the right hand side of equation 
(43) is written in a manner identical to that of the elliptic equations in 
equation (35). 

Stability of the Fourier approximation of the heat equation is easy to 
prove and is discussed in Gottlieb, Hussaini, and Orszag [6]. The more com- 
plicated case is equation (43) above. Kreiss and Oliger [17] have proved 
stability with two different treatments of the first order term. The first 
treats it in "skew symmetric form", that is, by writing 

Bu + Cu = -J- (Bu + (Bu) ) + (C - i B )u. (44) 

x 2 x x 2 x 


The second, discussed more fully in the next section, involves filtering the 
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first derivative to ensure stability, A convergence estimate using the skew 
symmetric form for the scalar equation with f = 0 is 

► 

T + 1 

Theorem 11: Let x 1, T > 0, and assume that u o € p v® , 2w ) • Then 

there exists a constant C such that the following estimate holds: 

j llu(t) - u c ll 0 _< C(1 + N 2 ) (1 _t)/2 {Hu 0 II t _ 1 + [Hu 0 ll 2 + (Hu 0 ll T + IIu 0 II t + 1 ) 2 ] 1/2 } (45) 

► 

for all t€[0,T]. 

I 

i Proof : See Mercier [10] Theorem 11.2. Here H p is defined in terms of 

distribution derivatives of periodic functions. 

) The convergence of Chebyshev approximations to parabolic equations on 

bounded domains has received a lot of attention recently. The spatial approx- 
imation for a Dirichlet problem will be exactly like that for the elliptic 
problem. Stability for the heat equation with non-constant coefficients was 
originally shown by Gottlieb [18]. Convergence estimates were worked out by 
Canuto and Quarteroni [19]. For the scalar heat equation 

u = a(x)u x € (-1,1) (46) 

t xx 

t with homogeneous boundary conditions u(-l,t) = u(l,t) =0 they show 

Theorem 12: Let a > 1/2 and S > a +2 and T > 0. If 
1 S 

u€L (0,T; H (-1,1)) then there is a constant C, Independent of N such 
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that 

Hu(t) - u (t)ll M < c N 2a+4-S 
c N — 


(47) 


where the norm II • IL, is the discrete norm derived from the Gauss-Lobatto 
N — — — — 

integration formula 


11 V " N 


N w.v(x )' 
y J J 

j-0 a( V 


(48) 


Neumann boundary conditions for parabolic problems can also be applied 
either explicitly or implicitly. For the implicit treatment, convergence is 
similar to that of the corresponding elliptic equation: 


Theorem 13: Suppose the solution to the differential equation (46) with 

2 in 

Neumann boundary conditions is regular to the extent that u€Lr (0,T;H™(-1 , 1 )) 

and the time derivative satisfies u^€L 2 (0,T;H in- *(-l,l)) for T > 0 and m 
■ t w 

> 5/2. U_ u 0 £H“(-l,l) then 

N— 1 « \ to 

llu(t) - « c <t)l 0>K + (u x - (u l)x ) Wj) 

< c i N2 '"f ll “o , »,» + et/2| / t,u "™ dt <49) 


Proof: See Canuto [16] Theorem 4.4 
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1.4.3 HYPERBOLIC EQUATIONS 

The study of the convergence of spectral approximations to hyperbolic 
equations is complicated by the fact that the straight-forward discretization 
of an equation of the form 


written as 


u t + a(x)u x + bu = 0 


3 u 3u 

7T + p N <a(x) nr> + bu c ' 0 


(50) 


(51) 


is often unstable. In this section, we will discuss the available theory of 
formally stable approximations. 

Fourier methods are stable if a(x) is of fixed sign. If a(x) in 
equation (50) is strictly positive and b = 0, the energy estimate for the 
approximation, (51) 


d 

Jt 


N-l u^x.) 
V c 3 

in ^7 


= 0 


(52) 


shows that the approximation is stable. If a(x) is zero at some point, how- 
ever, then this estimate is not valid and no general technique is available to 
show stability. 

Two basic approaches have been used to devise schemes which can be shown 
to be stable. The first, indicated in the last section, is to write the spa- 
tial derivative in skew-symmetric form. That is, instead of computing (51), 
one computes 


au c 

TT 


+ 1/2 P 


N^" 


a3u c 3(P N au) 

+ 51E i 


P nO + bU c 


0. 


(53) 
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Kreiss and Oliger [20], [17] showed that this discretization is stable. 
Mercier [10] examined the stability and convergence of the Fourier approxima- 
tion to the skew-symmetric equation 

u + v(x)u + (v(x)u) = 0 (54) 

u X X 

and showed that the error decay is spectral* 

Theorem 14: For t > 1 and T > 0, if the initial condition 

satisfies u(x,0)€H^ T (0, 2 tt ) then there is a constant C independent of N 
such that 

llu(t) - u (t)l < C(1 + N 2 ) (1 " x)/2 Uu n ll . (55) 

Proof: See Mercier [10], Theorem 9.1. 

Though approximations written in skew symmetric form are stable, there 
are objections to their use. The first objection is that they are less 
efficient since they have twice as many derivatives to evaluate. More 
important, conservation is lost when this is applied to conservation law equa- 
tions (such as the equations of gas-dynamics) for the computation of weak 
solutions. Tadmor [21] has examined the skew-self adjoint form of systems of 
non-linear conservation laws. They can be explicitly shown to be well-posed, 
but the conservation property is lost. 

The alternative to rewriting the equation in skew-symmetric form is to 
use the approximation of equation (51) and filter the solutions. Finite 
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difference solutions are often filtered by adding an explicit low order 
artificial viscosity. The goal of filtering Fourier spectral solutions is to 
do so without destroying the accuracy of the method. 

Two approaches for filtering Fourier approximations to guarantee 
stability have been suggested. The first was proposed by Majda, McDonough and 
Osher [12). In their method, the spectral derivative defined in equation (12) 
is modified by filtering the computed solution. For linear problems, this can 
be done efficiently by modifying the Fourier coefficients of the solution and 
using those new coefficients when the derivative is computed. Let 
p (x)ec°°(— it ,ir ) be a "filter function". Its values are zero near 9 = ± n 

and identically one in a neighborhood of 9=0. The Fourier coefficient 

* a 

u^ is replaced by p(2irk/N)u^ and this is used in equation (12) to 
compute the derivative. 

For smooth initial conditions, smoothing gives a stable approximation and 
spectral accuracy 


Theorem 15: 


u o ec 


the error satisfies the inequality 


llu(x,t) - u (x,t)l 
c 


< Ch' 
s — 


for all s,X 


depends on both 


Proof: See Majda, McDonough, and Osher [12], Corollary 1. 


If the solution is discontinuous, it is still possible to obtain spectral 
accuracy in the sense of equation (56) if the initial condition is properly 
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sraoothed. It is not enough to smooth the discrete Fourier coefficients of the 
initial condition with a filter whose support is enclosed within the support 
of p. Rather, it is necessary to use smoothed versions of the exact 
Fourier coefficients. 

A different approach to filtering was proposed by Kreiss and Oliger 
[17], Instead of filtering the solution with a predefined filter, they showed 
that linear stability could be obtained by smoothing the space derivative with 
a weak filter which depends on the smoothness of the function. They arbi- 
trarily split the frequency range of the solution into a high frequency range, 
|k| > Np and a low frequency range, |k| < Np The coefficients of the low 
frequency range are not modified at all. The coefficients of the high fre- 
quency range are modified only if they do not decay rapidly enough. Call 
v(x) = defined by equation (12) and define vj to be the derivative 
summing only the low frequency components |k| Np The modified 
coefficients of the derivative are defined to be w = Hv where 


w k = \ v k 


for | k | < 

if | k | > N 1 and | v R | < D II Vj R / | k | J (57) 

D ivl,/(|y| | k | ^ ) otherwise. 


Kreiss and Oliger prove the following stability theorem: 

Theorem 16: Suppose the coefficient a(x) in equation (50) is smooth so 
that its Fourier coefficients decay at a rate |k| The approximation 


H9u c 

u t + V a ■ 0 


3x 


( 58 ) 


-25- 


where the filter H Is defined by (57) is stable and converges if j = 8 > 2. 

Proof: See Kreiss and Oliger [17], Theorem 4.2. 

For linear problems, however, it is not clear that filtering is always 
needed. The fact that the energy method gives only a sufficient condition for 
stability means that equation (52) does not prove instability if a is not of 
one sign. For example, Gottlieb, Orszag, and Turkel [22] show stability in 
the usual sense of convergence as N °° of the scalar equation where a(x) 
= Asin(x) + Bcos(x) + C for arbitrary A, B, C. The numerical solutions do, 
however, grow in time - just as the exact solutions do. 

For non-linear problems, experience shows that filtering of the Fourier 
approximation is needed, particularly if there are discontinuities in the 
solutions. Hussaini, Kopriva, Salas, and Zang [51] discuss the application of 
these filtering methods and the choice of filters to a periodic transonic flow 
with a shock. 

Proofs of the stability and convergence of Chebyshev approximations have 
the added complication of the boundary conditions and the weight, w(x), which 
is unbounded at the endpoints. In particular, the case where a(x) changes 
sign makes it difficult to show stability. Gottlieb [18] has proved stability 
of the straightforward Chebyshev collocation for the special cases where a(x) 
= ±x . 

To show stability of Chebyshev approximations in general, the skew- 
symmetric form of the equations is needed. We will survey the convergence 
theory of Canuto and Quarteroni [23] for the special case of the hyperbolic 
boundary value problem with b(-l) > 0 
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u t + (bu) x + b Q u = f x€(-l,l), t€(0,T] 

u(-i,t) = o te(o,T) 

u(x,0) = Uq x€(-l,l). 

The further assumtion is added that 


(59) 


1/2 b + b - 1/2 bw w" 1 > 0 for x€(-l,l). (60) 

X O X — 

For the Chebyshev weight w(x) = (1 - x 2)-l/2^ th e use Q f integration by 
parts to get an energy estimate will give an unbounded boundary term evaluated 
at x = +1 (see Gottlieb and Orszag [!])• This has led to the use of a modi- 
fied weight and norm with which to prove stability and convergence. Let the 
new weight be w*(x) * (1 - x)w(x) so that w*(l) * 0, Then the following 
convergence estimate holds: 


Theorem 17: Suppose that u€L°°(0,T; H^(-l,l)) and b€ H®^(-1,1) 

w w 

for 8 > 2. Then the skew-symmetric Chebyshev approximation to (59) 
satisfies 


Hu - u II _ < CN 2 0 llull“ /ti e .. 

C L°(L 2 *) - L (H *) 

w w 


(61) 


Proof : See Canuto and Quarteroni [23], Theorem 2.3. Note: Their 

theorem actually allows for more general boundary conditions than we have men- 


tioned here 
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For computations which are not done in skew-symmetric form, such conver- 
gence estimates are not available in the general case. As indicated above, 
Gottlieb [18] has shown stability in some particular cases. However, Reyna 
[24] has shown that if b(x) is not strictly positive in the interval that a 
straight-forward Chebyshev approximation need not be stable. To stabilize the 
solutions he proposes the use of filtering. It is not sufficient, however, to 
simply smooth the Chebyshev coefficients. Rather, he shows that stability can 
be proved if Legendre coefficients are computed from the Chebyshev ones, the 
Legendre coefficients are smoothed, and then transformed back. 

1 The stability and convergence of Chebyshev approximations to the hyper- 

| bolic initial-boundary-value problem for systems 

i 

U t “ Au x -i<x£ i, t > 0 (62) 

i 

where u is an n-vector and A is a constant matrix has recently been proved 
by Gottlieb, Lustman, and Tadmor [25], [26]. Because this system is hyper- 
bolic, the matrix A can be assumed to have been diagonalized to 

TA 1 0 

A = 

.0 A 

where A* < 0, A** > 0 are diagonal matrices. 

Boundary conditions for which this system is well posed are of the form 

u X (-l,t) = Lu XI (-l,t) + g X (t) 
u n (l,t) = R u X (l,t) + g XX (t) 



( 63 ) 
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where and u^ represent the partition of u into inflow and outflow 

components (see Kreiss and Oliger [5]). Under the assumption that the bound- 
ary conditions are dissipative, the standard Chebyshev collocation is stable: 

Theorem 16: Under the assumption that |r| • |l| _< 1 - 6 < 1, the 

Chebyshev collocation method is stable for the system (62) with boundary con- 
ditions (63) in the sense that there exists a weighting pair w(x) and con- 
stants q and _> 0 such that for all s with Re s = n > Uq 

(n - n 0 )iiu c (x,s)n 2 < CN 2q |g(s)| 2 

where u_ and g are the Fourier transforms of u„ and g. 

Proof: See Gottlieb, Lustman, and Tadmor [25]. 


2. SOME APPLICATIONS OF SPECTRAL COLLOCATION METHODS 

In this section, some recent developments in the application of spectral 
methods to problems in fluid mechanics are surveyed. Much current emphasis 
has involved making spectral methods more efficient and more applicable to 
problems with complicated geometries. This has lead to the development of 
spectral multidomain methods which eliminate the need for global mappings and 
to the development of iterative techniques for the rapid inversion of the full 
matrices which occur when implicit time discretizations are used. 
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2.1. METHODS FOR ELLIPTIC AND PARABOLIC PROBLEMS 

2.1.1. SPECTRAL MULTIDOMAIN METHODS 

Spectral multidomain methods have been developed in order to avoid the 
need for global mappings required by spectral methods in problems with compli- 
cated geometries. A complicated domain can be subdivided into several subdo- 
mains and individual spectral discretizations can be applied to each subdo- 
main. For elliptic and parabolic problems, for handling the interfaces, early 
work considered explicit enforcement of continuity (e.g. Orszag [27] and 
Morchoisne [28]). More recently, spectral element discretizations and en- 
forcement of global flux balance have been used. The spectral element methods 
retain the accuracy of spectral methods in the context of a geometrically 
flexible finite element formulation. Global flux conservation has been used 
effectively when the mappings and/or domain sizes vary radically across inter- 
faces. 

Consider first the solution of the (second-order, self-adjoint, elliptic) 
Helmholtz equation, 

7 2 u - X 2 u = f in D (64) 

with Dirichlet boundary conditions on the domain boundary, 3D. Following 
the lead of finite element techniques, the spectral element algorithm [29, 30] 
proceeds by recognizing the equivalence of (64) to maximization of the follow- 
ing variational form, 


J d {-Vu»7u/ 2 - A 2 u 2 /2 - uf}dx. 


max 

ueH 1 


(65) 
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The variational form, (65), is preferable over the differential statement, 
(64), in that it requires less continuity of candidate solutions. 

The spectral element discretization proceeds by breaking up the computa- 
tional domain, D, into general quadrilateral elements. Within a given element 
k, the solution, geometry, and data are then expanded as tensor product 
Lagrangian interpolants through a set of specified collocation points. For 
instance, in two space dimensions, the solution u in element k is 
represented as, 

u k (r,s) = l l uf h . (r)h (s) (66a) 

-l -i 1 J 


h 

m 


(z 


n 


) 


(66b) 


where r and s are the local elemental coordinates, the h^ are the 
Lagrangian interpolants, the z n are the collocation points, and 6 is 

the Kronecker delta symbol. All summations run from 0 to N, where N is the 
order of the Lagrangian interpolants in each element. 

The expansions (66a) are then inserted into (65), and the functional 
rendered stationary with respect to arbitrary variations in the nodal values, 
u^ * Direct stiffness summation [31] (which recognizes that the global 
approximation space must be C®) is then used to assemble the elemental equa- 
tions into the system matrix. It should be noted that, as regards the treat- 
ment of elliptic and parabolic equations, the "spectral element" recipe pre- 
sented here is very similar to earlier "p-type finite element" methods [32] 


and the "global element" method [33]. 


It is clear from the above representation, (66), that the global inter- 
polant space is only C^, that is, that the approximation space suffers dis- 
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continuities in derivative at elemental boundaries. Although this may appear 
to violate the basic smoothness required of spectral methods, this is not the 
case due to the fact that the variational formulation, (65), is used rather 
than the (unintegrated) Galerkin weighted-residual form. In particular, in 
the absence of "variational crimes", the spectral element method can be shown 
to achieve exponential convergence to smooth solutions as N, the order of 
(fixed) elements, is increased. For nonsmooth solutions (e.g., corner-induced 
singularities) , high-order convergence is more difficult to obtain, however 
refinement techniques have been developed for the p-type finite element method 
[32]. 

Variational crimes take the form of numerical quadrature errors and in- 
terpolation of boundary data. (Nonconforming elements are not considered.) 
In order to insure that these errors do not dominate the approximation errors, 
it is important to correctly choose the collocation points of the Lagrangian 
interpolants . Earlier work on spectral elements used the Chebyshev colloca- 
tion points, as they are simple to evaluate and amenable to fast transform 
techniques. However, as the variational formulation (65) has essentially a 
unity weighting, it appears that a better choice is the Legendre-Lobatto 
points from the point of view of accuracy and efficiency of numerical 
quadratures [6, 34], Although Legendre polynomials are less convenient than 
Chebyshev polynomials, are subject to round-off errors for high-order expan- 
sions, and cannot be "fast transformed", none of these objections are 
particularly oppressive for the relatively low-order expansions used in 
spectral element methods. 

As an example of the accuracy of a Legendre spectral element [34] (see 
[29, 30] for extensive discussion of Chebyshev-based techniques), consider the 
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problem 

V 2 u =0 in D (67a) 

where D is the domain defined in Figure 1, x€[0, 1], ye[0, 1 + simrx]. 
Dirichlet boundary conditions are imposed such that the solution to the prob- 
lem is given by, 

u(x,y) = sin(x)e - y. (67b) 

The L error for the spectral element mesh shown in Figure 1 is plotted in 
00 

Figure 2. As expected from the analytic nature of the solution (67b) in the 
complex plane, exponential convergence is achieved as the order of the 
elements is increased. 

As another example of elliptic problems, consider the moving-boundary 
Stefan problem [34], given by 


V 2 0 =0 in D , D 2 (68a) 

0 = 1/2 on 8D (68b) 

V0* n|_ + 3V0*n| + = 0 on 9D ] . (68c) 

V0 • n = 0 on 3 D q (68d) 

0 = 1 + 1/2 cos2ttx on 3D^ (68e) 

0=0 on 3D 2 , (68f ) 
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where Dj, D 2 , 3Dj, 3 Dq, 3Dj, and are defined in Figure 3. Here the 

evaluation of _ and + refer to the Dj and D 2 sides of 3Dj, 
respectively. In point of fact, the time-dependent (parabolic) version of 
(68) was solved, approaching the steady-state only as t +• »; since the 
solution of parabolic equations involves at each time step the solution of an 
elliptic equation of the form (64), this aspect of the problem does not 
warrant separate discussion. 

Solution of the Stefan problem (68) illustrates several aspects of the 
spectral element method. First, since the interface 3D^ is unknown and 
general, it demonstrates the ability to handle complex geometry. Second, 

though the solution suffers a discontinuity at 3D^ the method has the 
ability to resolve certain non-homogeneities without losing "spectral 
accuracy". Figure 4 shows the interface position obtained with a Legendre 
spectral element method using a two-element mesh. In Figure 5, the associated 
temperature (0) distribution is given. High accuracy can be achieved with 
very few points. 

It is critical that the spectral element schemes not only be accurate, 
but also efficient as regards work required for a given level of accuracy. 
The key to the computational efficiency of the techniques is the sum factori- 
zation which follows from the tensor product representation, (66). For in- 
stance, a typical elemental term in a two-dimensional Chebyshev spectral 
element equation is of the form, 



(69) 


where h^, are defined as In (66), and all subscripts range from 0 to 
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N, the order of the polynomial space in each co-ordinate direction. Naive 
evaluation of this sura gives an operation count of O(N^), and O(N^) in 
three dimensions. This sura factorization is at the heart of both direct 
solvers using static condensation and fast eigenfunction solvers [35] and 
iterative solvers using conjugate gradient algorithms [36]. 

Another approach to handling domain interfaces was taken by Macaraeg and 
Streett [37], [38]. Within subdomains, the usual collocation procedure 

described in Section 1.2.2 is used. The interface values are computed by 
requiring that the solution be continuous and that the global flux be 
balanced. As an example of the procedure, consider the equation 

G(u) = F (u) - vu * S(u) 
x xx 

u(-l) = a (70) 

u(l) = b 


where an interface is placed at x - Integration of (70) from -1 to 
1 and the requirement that the jump in the flux [G] be zero at the interface 
yields 


G(u) . + / S(u)dx * G(u) - / S(u)dx. 
x=-l J . x=l J 

~1 x. 


(71) 


Numerical experiments show that spectral accuracy is retained. In two dimen- 
sions, the method has been used to solve Laplace^s equation with discontinuous 
boundary conditions. 
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2.1.2 ITERATIVE SPECTRAL METHODS 

For evolution problems, explicit time-stepping can be extremely ineffi- 
cient. This is because the typical time-step limitation for spectral methods 
is proportional to 1/N Z for the advection equation and 1/N 4 for the dif- 
fusion equation (where N is the number of modes) [39]. Hence, implicit time- 
stepping becomes a necessity. This results in a set of algebraic equations 
which are, in general, amenable to iterative solution techniques only. Also, 
elliptic equations governing practical problems virtually require implicit 
iterative techniques. Since the condition number of the relevant matrices are 
large, preconditioned iterative schemes including raultigrid procedures are the 
attractive choices. In this section, the fundamentals of iterative spectral 
methods are discussed with reference to an elementary example. 

For the purpose of illustration, consider the equation, 

u x = f, (72) 

periodic on [0,2ir]. For the Fourier method, the standard choice of colloca- 
tion points is given in Equation (11). 

The Fourier collocation discretization of the equation (72) may be 
written 

LU = F, (73) 

where U = (u Q , Uj, ..., u ), F = (f Q , f , f p, and L = C -1 DC. 

Here C is the discrete Fourier transform operator, C ^ the inverse trans- 
form, and D the diagonal matrix denoting the first derivative operator in the 
Fourier space. Specifically, 
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and 



-2irik 

e 


(j-N/2) 

N 

9 



i(j - N/2) 
- 0 


j,k = 0, 1, ..., N-l 

for j = 1,2, . . . , N-l 
for j = 0 


(74) 


The eigenvalues of L are X(p) = ip, p = -N/2 + 1, ••• , N/2 - 1, and the 
largest one grows as N/2. A preconditioned Richardson iterative procedure for 
solving Eq. (73) is 

V <- V + tolT 1 (F - LV) (75) 


where the preconditioning matrix, H, is a sparse, readily invertible approxi- 
mation to L. An obvious choice for H is a finite difference approximation 

to the first derivative. With the various possibilities for Lp D , the 
eigenvalue spectrum of L,I1 l is given in Table I. Apparently, the staggered 
grid leads to the most effective treatment of the first derivative. This kind 
of preconditioning was successfully used in the semi-implicit time-stepping 
algorithm for the Navier-Stokes equations discussed in section 2.2 on Navier- 
Stokes Algorithms. The eigenvalue trends of that complicated set of vector 
equations are surprisingly well predicted by this extremely simple scalar 
periodic problem. 

Next, consider the second order equation 

-u = f on [0,2 tt] (76) 

xx 


with periodic boundary conditions. A Fourier collocation discretization of 
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this equation is the same as Eq. (73) except for the diagonal matrix D which 
represents now the second derivative operator in the Fourier space. 

D j j = I" (j -f) 2 > J - 1, 2, N-l . (77) 

(O, j = 0 

O 

The eigenvalues of L are X(p) = p , p = -N/2 + 1, N/2-1. To make 

the case for the multigrid procedure (consisting of a fine-grid operator and a 
coarse-grid correction) as a preconditioner, assume H to be the identity 

matrix I in the iterative scheme (75). The iterative scheme is convergent if 
the eigenvalues, (1 - rnX), of the iteration matrix [I - wL] satisfy 


1 - uX| < 1. 


Each iteration damps the error component corresponding to X by a factor 
v(X) = 1 1-wX | . The optimal choice of X is that which balances damping of 
the lowest-frequency and the highest-frequency errors, i.e.. 


(1 - o»X ) = - (1 - ojX ) 
max min 


(78) 


This yields 


U SG Tx + X ) * 
max min 


(79) 


and the spectral radius 


^max ^min^ 

u sg = Tx + x ) * 

max min 


(80) 


In the present instance, X = N z /4, X 
F 9 max 


min = 1, and thus y gG « 1 - 8/N . 
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2 

This implies order N iterations are needed for convergence. This poor 
performance is due to balancing the damping of the lowest frequency eigenfunc- 
tion with the highest-f requency one. The raultigrid procedure exploits the 
fact that the lowest-f requency modes (|p| < N/4) can be damped efficiently on 
coarser grids, and settles for a relaxation parameter value which balances the 
damping of the mid-frequency mode (|p| = N/4) with the highest-f requency one 
( | p | * N/2). Table II provides the comparison of single-grid and multigrid 
damping factors for N=64. The high frequencies from 16 to 32 are damped 
effectively in the multigrid procedure, whereas the frequencies lower than 16 
are hardly damped at all. But then some of these low frequencies (from 8 to 
16) can be efficiently damped on the coarser grid with N=32. Further coarser 
grids can be employed until relaxation becomes so cheap that all the remaining 
modes can be damped. In concrete terms, the ingredients of a raultigrid tech- 
nique are a fine-grid operator, a relaxation scheme, a restriction operator 
which interpolates a function from the fine grid to the coarse grid, a coarse- 
grid operator, and a prolongation operator interpolating a function from the 
coarse grid to the fine grid. The fine grid problem for the present example 
may be written 

L f U f = F f . (81) 

Let denote the fine-grid approximation. After the high-frequency content 

of the error has been sufficiently damped, attention shifts to the 

coarse grid. The coarse-grid problem is 

L C U C = F c (82) 


where 
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F C = R [F f - L f V f ] , 

R being the restriction operator. After a satisfactory approxmation, V c is 
obtained; the coarse-grid correction (V c - RV^) is interpolated onto the 
fine grid by the prolongation operator P, yielding the corrected fine-grid 
solution 

V f <- V f + P (V C - RV f ) (83) 

i 

i 

The details of spectral multigrid techniques are furnished in [40]. Spectral 
multigrid techniques have been used to solve a variety of problems including 
the transonic full potential equation [41, 42]. Additional applications of 
spectral methods to compressible flows are described in [42]. 

i 

i 

| 2.1.3 Convection-Dominated Flows 

A model for convection-dominated flows is the viscous Burger^s equation, 

[ 9 

+ (u )^/2 = vu^ u(x,t=0) = -simrx, (84) 

with boundary conditions u(-l) = u(l) =0, and "small dif fusivity ," 
v = ,01/ir [43]. The solution to this problem develops a near shock. This 

near shock is characterized by the time at which the derivative at the origin 
attains a maximum value, t max , and the value of its maximum derivative, 

|8u/3x| . The convective term is clearly dominant for short times, however 

the diffusion terra insures that the solution will be smooth. This convection- 
diffusion balance is a good model for the kind of phenomena that arise in 
solution of the incompressible Navier-Stokes equations. The critical nuraer- 
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ical issues are numerical dispersion and diffusion. The former leads to in- 
correct propagation speeds of the shock affecting t max . The latter leads to 
smearing of the shock affecting l^ u ^ x l m ax # 

This problem has been solved by a variety of methods, including the spec- 
tral element method [43] and the explicit flux balancing method [37]. The 
spectral element calculations have used Crank-Nicolson in time on the diffu- 
sion term and the resulting Helmholtz equation in space was solved using the 
variational methods presented in Section 2.1.1. The convective term was 
handled by explicit third-order Adams-Bashf orth. Four elements were used 
covering the intervals [-1., -0.05], [-0.05, 0], [0., 0.05], [0.05, 1.] which 
cluster points around the location of high function variation. Macaraeg and 
Streett [37] used three subdomains with their flux conservation method. Table 
III presents a comparison of various methods on this model problem. 


2.2. INCOMPRESSIBLE NAV IER-S TORE S EQUATIONS 

This section is devoted to a description of algorithms for the solution 
of the incompressible Navier-Stokes equations in primitive variable form. The 
algorithms are based on methods discussed in the previous section in the 
simplest context. For example, simulation of instability and transition to 
turbulence in a flat-plate boundary layer have used iterative methods 
described in section 2.1.2. The spectral element method has been used for a 
variety of flow computations, including the problem of flow past a cylinder. 

The Navier-Stokes equations in the so-called rotation form are 

q^ - q x a) + V • (pVq) - VP 


in D 
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and 


V • q « 0 
q(x,0) = q Q (x) 

q = g 


in D (85) 

in D 

on 3D 


where q = (u,v,w) is the velocity vector, u = V x q the vorticity, 

P - p + 1/2 | q | the total pressure, y the variable viscosity, D the in- 
terior of the domain, and 3D its boundary. In the stability and transition 
problems, the domain D is cartesian and serai-infinite: periodic in the two 
horizontal directions (x,z), and bounded by a wall at y=0. Fourier colloca- 
tion can be used in the periodic directions (x,z) and Chebyshev collocation is 
used in the vertical (y) direction. The collocation points in the periodic 
directions are given by a relation similar to Eq. (11). The vertical extent 
of the domain 0 < y < 00 is mapped onto -1 < g < +1. The velocities are 
defined and the momentum equations enforced at the points 


= cos (|i), j = o, 1, N y (86) 

The pressure is defined at the half points 

5 1 = cos [ * ] , j - o, 1 N -1 (87) 

j+j y 

and the continuity equation is enforced at these points. The staggered grid 
avoids artificial pressure boundary conditions, and precludes spurious pres- 


sure modes 
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After a Fourier transform in x and z, the temporal discretization 
(backward Euler for pressure. Crank— Nicolson for normal diffusion, and third 
or fourth-order Runge-Kutta for the remaining terms) of Eq. (85) leads to 

[I - MDM] Q + At A q VII = Q c (88) 

- A + V • Q = 0 

where 


Q 



~n+l 

<1 




~n+l 

% 


} 


n 


r ~n+l ^n+1 ~n+l x 

ip l/2 * p 3/2 p N-l/2 > 


V 



ik x> 


(89) 


M is the Chebyshev derivative operator, D the diagonal matrix with l/2yAt 
as its elements, and Aq is the interpolation operator from the half points 
to cell faces, A + vice versa. Obviously, the equations for each pair of 
horizontal wave number (k^ , k^) are indepenent and they can be written as 
the system 


LX = F 


(90) 


where X = [Q, II]. The iterative solution of this equation is carried out by 
preconditioning the system with a finite difference approximation on the 
Chebyshev grid, and applying a standard iterative technique such as 
Richardson, minimum residual or multigrid [44]. 
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The method described above solves the implicit equations together as a 
set. The extension of this method to the more general cases of interest such 
as those involving two or more inhomogeneous directions is not straightfor- 
ward. An alternative is the operator-splitting technique or the fractional 
step scheme [45]. This method yields implicit matrices which are positive 
definite and are easily amenable to iterative methods. In the first step, one 
solves the advection-dif fusion equation 


q t = q x a) + V • (yVq ) 


(91) 


subject to the initial and boundary conditions 


q (x»t n ) = q(x,t n ). 


* * 

q = g on 3D. 


(92) 


Note that g has yet to be defined. In the second step, one solves for the 
pressure correction 

•A.4. 4»4. 

(93) 

(94) 

subject to the conditions 


** ** 
q t = -vp 

** 

V • q =0 


**. * . * * 
q (x,t ) = q (x,t ) 

** 

g • n = g • n 


in D 
in 


(95) 


where n is the unit normal to the boundary. 


Further, the tangential 
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component of the Eq. (95) holds on the boundary, i.e., 

+ * A JLJU A 

q t • x = - VP • t in 3D (96) 

where t is a unit tangent vector to the boundary* Now g* is defined 
[45] as (using Taylor expansion in t) 

g* • n = (g n + At g£) • n (97) 

g • x = [g n + At (g* + VP n )] • x . 

Eq. (91) is discretized in the usual spectral collocation manner. After a 
temporal and spatial discretization of Eq. (93), the boundary conditions are 
built into the relevant matrix operators, and then a discrete divergence is 
taken. This results in a discrete Poisson equation (with as many algebraic 
equations as unknowns) for pressure, which can be solved by standard iterative 
techniques including the raultigrid method. 

Spectral element methods have been applied to the incompressible Navier- 
Stokes equations (85). In addition to (98), the uncoupled (passive) or 
coupled (natural convection) energy equation is also often of interest. The 
time discretization used for the Navier-Stokes equations is either a Green's 
function technique [29] or an operator splitting scheme [30]. Both of these 
methods reduce (85) at each time step to an initial convective step, followed 
by a Stokes problem consisting of a sequence of Poisson and Helmholtz equa- 
tions. The spatial discretizations discussed above in Section 2.1.1 are 
directly applicable to these Navier-Stokes subproblems. 
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Spectral element methods have been applied to the simulation of numerous 
flows in the Reynolds number range 0 < R < 1500 [36, 46-47], An example is 
provided by the classical problem of flow past a cylinder. Results are pre- 
sented here for R = 100, based on freestream velocity and cylinder diameter, 
for times sufficiently large that the flow has reached a steady-periodic 
state. Figure 6 shows the spectral element mesh used, and Figures 7 and 8 
show the streamlines and isotherms, respectively, at one time in the periodic 

flow cycle. The thermal boundary conditions are T = T far from the 

00 

cylinder, and T = T w on the cylinder surface. The isotherm pattern clearly 
reveals the spatial structure of the von Karman vortex street. Note the mini- 
mal numerical dispersion in the scheme, as evidenced by the clear identity of 
the shed packets of fluid and the absence of trailing waves in the wake. More 
details of these cylinder calculations, as well as comparisons with previous 
numerical work and experiment, can be found in [36], 


2 . 3 HYPERBOLIC EQUATIONS 

Here, the application of spectral methods to the solution of inviscid 
compressible flow problems is surveyed. Methods for such problems are not 
nearly so advanced as those for incompressible flows. The survey is limited 
to methods for the solution of the Euler equations of gas-dynamics governing 
some flows of aerodynamic interest. For the solution of the full potential 
equation for transonic flows, see Streett, et al. [42]. 

The Euler equations of gas-dynamics are a coupled system of nonlinear 
hyperbolic equations which (in one dimension) are usually written in the con- 


servative form 
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u t + F x (u) = H * (98) 

Typically, spectral discretization in space and explicit finite difference 
discretization in time are used. The discontinuous solutions of this set of 
equations have been obtained in the case of a shock tube (Gottlieb, Lustraan, 
and Orszag [48], Cornille [49]), quasi-one-dimensional flow in a nozzle (Zang 
and Hussaini [50]) and for the astrophysical problem of shocked flow in a 
galaxy (Hussaini, et al. [51]). 

The astrophysical problem is the most challenging one-dimensional com- 
pressible flow problem for which shock capturing has been attempted with a 
Fourier spectral method. It contains a strong shock and an adjacent strong 
expansion. Unlike problems with weak shocks and expansions, it was necessary 
to apply strong filtering to stabilize the numerical solution. The result of 
this drastic filtering was a reduction of the order of accuracy. Even in the 
smooth parts of the solution away from the shock, the accuracy was only first 
order. In view of the extra work involved to compute the spectral approxima- 
tions, it is not clear that spectral methods with filtering are a viable 
alternative to finite difference methods when strong shocks are captured. 

An alternative to capturing shocks is to treat them as boundaries. In 
this case, it is possible to compute the solutions using the nonconservative 
form 

Q t + AQ X = E (99) 

along with an ordinary differential equation for the motion of the shock. A 
number of two dimensional shock-fitted solutions are described in Hussaini, et 
al. [52]. These solutions include a shock/ turbulence interaction. 
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Shock/ vortex interaction and supersonic flow past a cylinder. When shocks are 
fitted, spectral methods do indeed outperform typical second order finite 
difference methods, as long as the solution is adequately resolved. Kopriva, 
et al. [53] compared the performance between MacCormack's method and the spec- 
tral collocation method for the shock/plane wave interaction problem and for 
the Ringleb problem. A comparison of the accuracy of the finite difference 
method vs. the spectral method is shown in Table IV. 

A multidomain method for the nonconservative form of the Euler equations 
suitable for use with shock-fitting has been described by Kopriva [54]. In 
each subdomain, the usual collocation method (Section 1.2.2) is applied. At 
Interfaces, however, a weighted average of the derivatives is used. In one 
dimension, 

Q t + aLQ k + A \ R = E (100) 

where denotes the solution vector at the interface and the derivatives 
superscripted with the L and R denote the left and right computed spectral 
approximations. For consistency, A^ + - A. The weighting corresponds to 
an upwind approximation 

A L = 1/2(A + |a| ) A R = 1/2(A - |a| ) 

where |a| = Z | A | Z — 1 and Z is the matrix of right eigenvectors. For many 
applications, this can be simplified by replacing |a| by a diagonal approxi- 
mation | A | « X I where M min < X < M raax is an approximation to the 


eigenvalues of A. 
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Subdivtdtng the domain retains spectral accuracy. Table V shows the 
performance of a two domain computation of a one dimensional 2x2 

hyperbolic system 


u 


1 2" 


"u~ 


+ 




- V r- 


.2 1. 


V 


( 102 ) 


from [54] with solution u and v for an equal number of points on each side 
of the interface. 

For a given number of grid points, it is possible to obtain solutions 
with a multidomain spectral method which are significantly better than the 
single domain method. For the two dimensional nonlinear Ringleb problem 
computed by Kopriva [54], Table VI shows the effect on the error of a four 
domain division in which the position of the streamwise interface is varied. 
The accuracy is best when rapid changes in the solution are best resolved. 

The shock-vortex interaction problem described by Kopriva [55] provides 
another example of the advantage of a multidomain method over a single domain 
method. A two dimensional region between the shock and an arbitrary upstream 
boundary is mapped onto a square. The shock moves downstream where it en- 
counters a vortex. The interaction of the shock and the vortex creates a 
circular sound wave centered on the vortex. Because the physical domain is 
continually increasing in size, the resolution of the solution decreases with 
time. 

The single domain solution to the shock-vortex problem cannot be computed 
without added smoothing. Figure 9a shows the pressure contours with no 
smoothing. The numerical oscillations in the pressure are of the same order 
as the sound pressure wave created by the interaction. If the region between 
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the shock and the upstream boundary is subdivided and some of the subdomains 
are allowed to move with the shock, smoothing is not required. Figure 9b 
shows the pressure contours of a two domain calculation with the same number 
of grid points in the horizontal direction. The horizontal numerical oscilla- 
tions are no longer present and the sound pressure wave is clearly visible. 
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Figure 1 . 


Figure 2. 


Figure 3. 


Figure 4. 


Figure 5. 


Figure 6. 


Figure 7. 


Figure 8. 


Figure 9 


FIGURE CAPTIONS 


The domain and Legendre spectral element discretization used for 
solution of the Laplace equation described in the text. Although 
results are given for a single isoparametric element, similar re- 
sults are obtained with multiple elements. (Legendre results due 
to E. M. Ronquist.) 


A plot of the L error as a function of the number of points 
in one direction for solution of Laplace's equation in the domain 
shown in Figure 1 . 


Description of the domain and boundaries for the Stefan problem 
presented in the text. (Stefan problem results due to E. M. 
Ronquist. ) 


Spectral element prediction for the position of the interface, 

SDj, for the Stefan problem described in the text. The spectral 
element mesh uses two elements, one each in each phase. 


Temperature (0 ) distribution for the Stefan problem described 
in the text. Note the discontinuity of slope at the interface, 


Spectral element mesh used for simulation of flow past a cylin- 
der. Note the flexible resolution afforded by the elemental 
decomposition. (Flow past a cylinder results due to G. E. 
Karniadakis.) 


Instantaneous streamlines of the cylinder flow at a Reynolds 
number of R = 100. 


Instantaneous isotherms of the cylinder flow at a Reynolds number 
of R = 100 (Prandtl number of unity). The von Karman vortex 
street can be clearly seen in the temperature distribution behind 
the cylinder. 


Pressure contours for a shock/ vortex interaction. (a) Single 
domain calculation, (b) Three vertical domain calculation. 
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Figure 9a 
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Table 1. Preconditioned Eigenvalues for One— dimensional First Derivative 

Model Problem 


Preconditioning 


Eigenvalues 


Central Differences 


kAx 

sin(kAx) 


High Mode Cutoff 


kAx 

sin(kAx) 


IkAxI < (2n/3) 


(2ir/3) < kAx < ir 


One-sided Differences 


-i(kAx/2) kAx/2 

e sin((kAx)/2) 


Staggered Grid 


(kAx) /2 
Ln( (kAx) / 


m 

Single-Grid 

Multigrid 

1 

.9980 

.9984 

2 

.9922 

.9938 

4 

.9688 

.9750 

8 

.8751 

.9000 

12 

.7190 

.7750 

16 

.5005 

.6000 

20 

.2195 

.3750 

24 

.1239 

.1000 

28 

.5298 

.2250 

32 

.9980 

.6000 
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Table III. Comparison of Methods for Solution of 
Burger's Equation (from Ref. 43, 37) 


Method 

Interval 

Hg 

IT • t 

max 

N/M 

At»TT 

- Fourier Galerkin 

[-1,1] 

151.942 

1.6035 

682/1024 

5.10" 4 



142.665 

1.60 

682/1024 

10 -2 



148.975 

1.603 

170/256 

5.10 -4 



142.313 

1.60 

170/256 

10 -2 

- Fourier Pseudospectral 

[-1,1] 

142.606 

1.60 

256/256 

10 -2 



144.237 

1.60 

128/128 

10" 2 

- ABCN collocation 

[-1,1] 

145.877 

1.60 

512 

5.10“ 3 

+coordinate transform 

[-1,1] 

152.123 

1.60 

64 

10 -2 

- Spectral Element 

[-1,1] 

152.04 

1.6033 

16 x 4 

10“ 2 /6 

- FD 

[-1,1] 

150.1 

1.63 

81 

10" 2 

- Chebyshev 






ABCN spectral 

[0,1] 

152.05 

1.60 

64 

1/300 

Rosenbrook spectral 

[0,1] 

151.998 

1.60 

64 

10“ 2 


[0,1] 

150.144 

1.60 

32 

10" 2 

ABCN collocation 

[0,1] 

152.126 

1.60 

64 

10~ 2 

- Flux balance 

[-1,1] 

152.00011 




- Analytical 


152.00516 

1.6037 
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Table IV. Maxima Error in p for MacCormack and Spectral Computation 

of Transonic Ringleb Flow 


Grid 

MacCo rmack 

Spectral 

9x5 

2.6 x 10” 2 

2.2 x 10" 2 

17 x 9 

i.i x ice 2 

1.9 x 10 -3 

33 x 17 

3.2 X 10“ 3 

5.0 x 10 -5 

Table V. Solutions 

to (102) with Equal Number of 
of the Interface 

Point 8 on Each Side 

N 

Error in u 

Error in v 

8' 

1.57 x 10 -2 

1.49 x 10 -2 

16 

4.15 x 10 -6 

4.86 x 10 -6 

32 

1.91 x 10~ 9 

1.91 x 10 -9 


Table VI. Effect of Streanri.se Mesh Distribution on Ringleb Calculation 


Grid 

Division 

Maximum Error 

8 + 8 

0.45 + 0.55 

7.8 x 10~ 4 

8 + 8 

0.50 + 0.50 

9.3 x 10 -4 

16 (SD) 

- 

1.9 x 10" 3 

10 + 6 

0.34 + 0.66 

1.2 x 10 -2 
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